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Abstract 

The statistics of earthquakes in a heterogeneous fault zone is studied analyti- 
cally and numerically in a mean field version of a model for a segmented fault 
system in a three-dimensional elastic solid The studies focus on the 

interplay between the roles of disorder, dynamical effects, and driving mech- 
anisms. A two-parameter phase diagram is found, spanned by the amplitude 
of dynamical weakening (or "overshoot" ) effects e and the normal distance L 
of the driving forces from the fault. In general, small e and small L are found 
to produce Gutenberg-Richter type power law statistics with an exponential 
cutoff, while large e and large L lead to a distribution of small events com- 
bined with characteristic system-size events. In a certain parameter regime 
the behavior is bistable, with transitions back and forth from one phase to the 
other on time scales determined by the fault size and other model parameters. 

The implications for realistic earthquake statistics are discussed. 
PACS numbers: 91.30.Px, 05.40.+j, 62.20.Mk, 68.35. Rh 
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I. INTRODUCTION 



The statistics of earthquakes has been a subject of research for a long time. One spec- 
tacular feature is the wide range of observed earthquake sizes, spanning over ten decades in 
earthquake moment magnitude (which is defined to scale as the logarithm of the integral of 
slip along the fault during the earthquake ||). Gutenberg and Richter p| found in the 50's 
that the size distribution of regional earthquakes follows a power law over the entire range of 
observed events. The exponent b of the power law distribution appears to be universal, i.e. 
it is approximately the same (within statistical errors and possible secondary dependency 
on the tectonic domain) for all studied regions. This type of power law distribution is called 
the "Gutenberg Richter" distribution. Recently, enough data has been collected to extract 
statistics on individual systems of earthquake faults, or more precisely on systems of narrow 
fault zones. Interestingly, it was found that the distribution of earthquake magnitudes may 
vary substantially from one fault system to another. In particular, Wesnousky and coworkers 
0] found that fault systems with highly irregular geometry, such as the San Jacinto fault 
zone in California, which have many offsets and branches, display "power law" statistics 
over the whole range of observed magnitudes. Not all fault systems, however, display a 
power law distribution on all scales up to the largest earthquakes. The available data Q 
indicate that fault systems with more regular geometry (presumably generated progressively 
with increasing cumulative slip) such as the San Andreas fault in California display power 
law distributions only for small events, which occur in the time intervals between roughly 
quasi-periodic earthquakes of a much larger "characteristic" size which rupture the entire 
fault. There are practically no observed earthquakes of intermediate magnitudes on such 
geometrically regular fault systems. Distributions of this type are called the "characteristic 
earthquake" distribution. 

In previous work it was demonstrated that a class of simple models of ruptures 
along a heterogeneous fault zone displays both types of behavior. The universal power law 
scaling behavior of earthquake statistics was seen to be due to an underlying critical point, 
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which becomes mean-field-like for fault geometries with more than two spatial dimensions. 
In the limit of weak dynamical effects, the mean-field approximation to the 2 dimensional 
fault provides a more appropriate approximation than, for example, traditionally studied one 
dimensional approximations to the models. In fact, exact results for the scaling exponents 
(up to logarithmic corrections) could be obtained from mean-field theory. The reason is that 
the elastic stresses along the fault are effectively long range (decaying like the inverse cube 
of the distance), such that in two and higher dimensions the fluctuations due to interaction 
with other points on the fault decrease as the fault size is increased - on long length scales the 
behavior becomes the same as that of a system with infinite ranged elastic interactions (up 
to logarithmic corrections in two dimensions). In other words, the upper critical dimension 
is equal to the physical dimension of the fault, which is 2 [§|§. (Some of the static mean-field 
exponents turned out to be the same as in other quasi-static models ||.) In the presence 
of small but nonzero weakening effects of amplitude e a critical rupture size (slipping area) 
n cr for "runaway" or "characteristic fault size" events was calculated perturbatively |J and 
was found to scale as 1/e 2 . Faults of larger area than this size are expected to display the 
characteristic earthquake distribution, with small events up to size n cr , and no events of 
intermediate size between n cr and the characteristic fault size events. For faults of smaller 
total area than n cr only the power law scaling region of the small events is seen, so the 
distribution is of the Gutenberg Richter type. 

In this paper we examine a mean-field model with a range of dynamical weakening effects 
from weak to strong, and different levels of disorder in the brittle properties. Specifically, we 
study the model of Ben-Zion and Rice [|IJ, which involves simple approximations of dynamic 
frictional weakening (similar to static versus dynamic friction), but replace the physical long 
range elastic interactions with infinite range interactions. In addition to exhibiting both 
"power law" and "characteristic" scaling of event sizes, this model exhibits the possibility 
of coexistence of these two types of behavior. That is, for a given set of model parameters, 
the system has two distinct persistent stationary states. In an infinitely large system it will 
depend on the initial conditions whether the system displays Gutenberg Richter or charac- 
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teristic earthquake type behavior. Faults of finite size can spontaneously switch from one 
state to the other on time scales that are exponentially large in system size. The switch- 
ing times (or "persistence times") are determined by nucleation processes from one state 
to the other, similar to flips back and forth at coexistence in finite thermally equilibrated 
systems. Many of the qualitative features seem to be sufficiently robust to be applicable to 
real fault zones. Interesting to note, such "switching" behavior appears to characterize long 
paleoseismic records observed along the Dead Sea transform fault system in Israel ||, and 
is compatible with other paleoseismic [10| and geologic JTTJ data. In addition, qualitatively 
similar switching has been recently found in regional models of disordered fault systems |12| . 



The remainder of this paper is organized as follows: In Section |I|, we define the model and 



provide a summary of the main results. In Section [HT| , we present a detailed analysis of the 
model along with comparisons with numerical simulations. In Section [TV], we compare our 
results with earlier studies of similar models and discuss their potential relevance to natural 
fault systems modeled as a narrow fault zone in a three dimensional elastic surrounding 
media. 



II. THE MODEL AND SUMMARY OF RESULTS 

Ben-Zion and Rice |] suggested that a heterogeneous fault system with offsets and 
branches may be represented by an array of discrete cells in a two dimensional plane, with 
spatially varying "macroscopic" constitutive parameters that model the heterogeneity of 
the original fault system. This model fault on the (x, z) plane can be considered as a 
collection of brittle patches mapped onto the interface between two tectonic blocks, which 
move with (small) relative transverse velocity t>x far away from the fault. In the simple 
realizations used in refs. 0>@], and here, (as in related models 0), the fault plane is segmented 
into N geometrically equal cells. In the mean field approximation of infinite range elastic 
interactions, the local stress on cell i is given by 
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Ti = J/Nj2(uj - Ui) + K L (vt - Ui) 
j 

= Ju + K L vt- (K L + J)u h (I) 

where Ui is the total fault offset of cell i in the horizontal (x) direction, u = Uj)/N, J/N 
is the elastic coupling between cells in the mean-field approximation, and Kl is the effective 
loading stiffness of the bulk material surrounding the fault patch. 

Initially, the fault is in a relaxed configuration, i.e. all stresses are less than a local static 
failure threshold stress t s ^. In the absence of brittle failures the stresses at the cells increase 
uniformly due to the external loading and t$ = KlV. As long as no cell reaches its failure 
threshold, Ui = everywhere. When the stress at a cell becomes larger than t 8 j, the cell 
slips by an amount 5ui = (r si — T a>i )/(K L + J), to reduce its stress from r S) i to an arrest 
stress r aj j. (The nonuniformity of failure and arrest stresses across the fault plane models 
the spatial heterogeneity of real fault zones [|IJ.) Consequently, during failure cell stresses 
change by [cf. Eq.(]l|) ] 

Sn = r a>i - T Sti , (2a) 
<Jr i = (c/iV)(r s , f -r a , i ), j ^ i, (2b) 

where c = J / {Kl + J) is a "conservation parameter" giving the fraction of the stress drop 
of the failing cell retained in the system after the slip. As pointed out in Refs. for 
fault zones with characteristic linear dimensions of O(L), the "loading spring constant" is 
Kl ~ 1/L, provided that the stress loading of the fault is either due to uniformly moving 
(creeping) boundaries or applied forces at distances of O(L) away from the fault plane. For 
the case iV = L 2 , (1 — c) ~ O (l/ v^a) • A value c < f for a large system would be physically 
realized if the external drive is closer to the fault than its linear extent. 

During the failure process, the slipped cell is assumed to be weakened by the rupture, 
such that its failure strength is reduced to a dynamical value r di = r Sj j — e(r s i — r aji ), 
with < e < I parameterizing the relative importance of dynamical weakening effects in the 
system. If the failure stress transfer brings other cells to their failure threshold, an avalanche 
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of cell failures, i.e., "rupture propagation", occurs according to Eqs.(0) until all cells are at 



stresses 7$ < r s i ||13|| . It is assumed that these avalanches happen on time scales short 
compared to the external loading time (i.e. v is adiabatically small), so that the external 
load is kept constant during an earthquake. In time intervals between earthquakes, all cells 
are assumed to heal completely, thus failure thresholds are reset to their static value r s>i and 
the external loading resumes until the next cell failure. 

In order to simplify notation, it is useful to introduce rescaled stress variables 



*a,i = 1 



Sd,i 







T a,i) 



(4) 



T,l ' ! 1 e(l- .s,„). (5) 



" T a ,i) 

such that cell failure always occurs when Sj = 1, and (s a ) = 0. (Here, () symbolizes 
averaging over all cells in the fault zone.) The arrest stress s a ^ is uncorrelated from cell 
to cell, and is picked once for each segment from a probability distribution p(s a ) with 
mean and a compact support (— W/2, W/2) of width < W < 2, which characterizes the 
heterogeneity of the fault system. [In our simulations we have used the parabolic distribution 
p(s a ) = 3(W 2 - As 2 a )/(2W 3 ), for -W/2 < s a < W/2 and otherwise.] Unless stated 
otherwise, the focus is on the small disorder limit W <C 1 and moderate values for e, which 
are considered fixed, and the properties of the system are analyzed as a function of varying 
conservation parameter c and system size N. (In the last section of the paper we discuss 
the effects of larger values of W as well.) The size of an earthquake refers to the number of 
cells that failed, (i.e. the "area" on the fault that slips in an earthquake). 

For iV — >• oo, depending on relative values of the system parameters, there are in general 
two possible steady-state distributions of cell stresses and of earthquake magnitudes. We 
refer to these as as "phases" : 

(A) The "Gutenberg-Richter" (G-R) Phase: This phase, possible in both regions 1 and 
2 of Fig. [TJ, is characterized by a distribution of earthquake sizes (n) of power law form. 
In infinite systems (N — > oo), it is given by 
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p^(n)«A/n 3/2 exp(— n/n c /), n < iV, (6) 

with a characteristic cutoff size n c / ~ 2(1 — c)~ 2 that diverges as c /* 1. (Finite size 
corrections are given in equation fll5|) below) . The stress s, at a given cell is independent of 
all others and is equally likely to take any allowable value, i.e., 

ds 

Prob(s < Si < s + ds) = , s a>i < s < 1. (7) 

Thus, the stress distribution in the G-R phase is given by 

P if \s) = r d Sa ^, (8) 

J —oo 1 — S a 

which is uniform and equal to p = ((1 — Sa)^ 1 ) in the interval (W/2) < s < 1 [see Fig. 0(a), 
inset]. In this phase each cell fails at most once during an earthquake, and therefore dy- 
namical effects are largely irrelevant. An infinitely large system which started in the G-R 
phase will remain it this phase. In finite systems N < oo, with parameters in region 2 of 
figure [I], however, a very large earthquake of size (1 — e)N/c or greater occasionally triggers 
dynamical effects that lead to a catastrophic "runaway" event in which all cells eventually 
fail and cause a substantial change in the stress distribution and subsequent evolution of the 
system, as outlined next. 

(B) The "Runaway" Phase: This phase is characterized by a quasi-periodic occurrence 
of system wide earthquakes in which all cells fail. As a result of dynamical effects, the stress 
Si in a cell immediately after such a "runaway" event is independent of other cells and is 
equally likely to take any value between its arrest stress and dynamical failure stress, i.e., 

ds 

Prob(s < Si < s + ds) = , s a ^ < s { < s dji . (9) 

The stress distribution is thus given by 

which is uniform and equal to p/(l — e) in the interval (W/2) < s < 1 — e — (eW/2) [see 
Fig. 0(b), inset]. The runaway event is followed by a quiescent period during which stresses 
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on the cells build back up from their dynamic failure value to near their static failure value. 
Subsequent small events are followed by the next runaway event, at which point the stress 
distribution is reset to Eg. flUf) . These background small events have a size distribution 
similar to events in the G-R phase, but with a different cutoff size. In an infinite system 
(for finite size corrections see equation ( |T8"D below): 



P, 



(r) 



(n) « A r n~ 3/2 exp(-n/n cr ), n < N, (11) 

n * 2(1 " £)2 (12) 
ncr (1-e-c) 2 ' 1 j 

which diverges as c \ (1 — e). However, this divergence is never observed, as the runaway 
phase becomes unstable against breakup into the G-R phase when c < c* = (1 + e) _1 
for the following reason: If the background small events during a cycle involve at least a 
fraction r c of the cells, the subsequent large event is unable to cause all of the cells to fail, 
since the cells that failed during background activity are farther away from their failure 
stress. This typically causes a spontaneous breakup of the bunched stress distribution and a 
resumption of the G-R phase. The fraction of cells needed to cause this breakup is given by 
r c = 1+e—C 1 = (c*) -1 — c~ l as is derived in section [TTI|. When c \ c*, the size of background 
events necessary to cause breakup vanishes and the runaway phase becomes unstable, i.e. 
for c < c* the G-R phase is the only persistent phase, regardless of the initial conditions. 
For c > c*, the G-R phase and the Runaway phase are both persistent in an infinite system. 
In an infinite system the initial conditions determine which one the system displays. In a 
finite system, however, exponentially rare earthquakes can lead to nucleation from the G-R 
phase into the Runaway phase and vice versa. Equations (|T6| ) and ( ^0|) below give estimates 



of the times spent in the respective phases between such nucleation (or "switching" ) events. 

III. ANALYSIS OF THE MODEL 

The results quoted above have been obtained by mapping earthquakes in the model to 
corresponding events in a stochastic process, which is approximated by a series of Bernoulli 
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trials [II] in order to be able to obtain analytical estimates for the various quantities of 
interest, such as distributions of earthquake sizes and persistence times for the two phases. 



A. Gutenberg-Richter Phase 

Let us first focus on the G-R phase. At some instant t immediately preceding a cell 
failure, consider the sequence jx„ = 1 — Si( n+ i)| , where i{n) is the index of the cell that 
has the nth largest stress in the system [See Fig. |3|]. For a large system, the stress gaps 
{Ss n = X n _i — X n } are (almost) independent of each other, drawn from an exponential 
probability distribution, i.e., Prob(£s n = s) = exp(-pNs), with p = ((1 — s a ) _1 ). For 
n ^> 1, X n resembles a biased random walk with a mean Hx{n) = n/(pN) and variance 



n/(pN) 2 . As long as dynamical effects are absent, the stress redistributed to each 



cell following the failure of the first n cells is given by a random variable Y n with mean 
fiy( n ) = nc/(pN) and variance <7y(n) ~ n(cW/pN) 2 <C cr\{n). A triggered earthquake 
can sustain itself only if the redistributed stresses exceed the stress gaps. Therefore, Z n = 
X n — Y n < during an earthquake and it immediately follows that the distribution of 
earthquake sizes for N ^> n ^> 1 are given in terms of the distribution of first passage times 
of biased random walks. Approximating the continuous probability distribution of the step 
sizes of {Z n } with a Bernoulli process (where steps of equal size are taken up or down with 



probability p and 1 — p, respectively), we can utilize results available for Bernoulli trials |14 



Prob (Zi < 0, < % < n; Z n = 0) = Prob(Z " 21, (13) 

n 

i.e. the probability for the first return to the origin after n steps equals the total probability 
of reaching the origin after n steps divided by n. Prob[Z n = 0] can easily be calculated |14[ 



for iV 3> n 3> 1. One obtains Eq.(D) with n c f = 2(u 2 + <7 2 )/u 2 , where u and a 2 are the 
mean and variance of the step size for the process Z n , respectively. Substituting the values 
/i = (1 — c)/(pN) and a 2 ~ l/(pN) 2 , the cutoff length is given by 

n cf » 2(1 + (1 - c)- 2 ) = — {l + O [(1 - c) 2 ] } , (14) 



where the last approximation is justified since treating Z n as a Bernoulli process is expected 
to yield relative errors of 0(/j, 2 /a 2 ). 

For finite-sized systems, when the fraction of failed cells r = n/N is no longer small, 
Eq.(H) needs to be modified since the stress gaps are not entirely independent: In order to 
correctly reflect the fact that = 1 — c to within 0(1/N), the Bernoulli process should be 
constrained to return to its mean value after N steps. This can be achieved by calculating 
the corresponding conditional probabilities: 

(/) _ Prob(Z„ = 0) Prob(ZAr_„ = 1 - c) 



n Prob(Z7v = 1 — c) 

k c::p f < l + 

n 3 / 2 \ n c f 



37^ ex Pi ~ h ( 15 ) 



which reduces to Eq.(|6|) in the limit n <C N. ( Af is a constant fixed by normalization.) 
Figure 0(a) shows the distribution of event sizes for numerical simulations of the model with 
N = 400, for values c = 0.6,0.7 and 0.8. (In all presented simulation results, W = 2/19 
and e = 0.5.) The continuous lines are one-parameter fits to the form (|T5|). The discrepancy 
between the fitted and theoretical [from Eq.fJTJ])] values of n c f is consistent with the expected 



relative error. 

As mentioned earlier, the failure of all the remaining cells becomes very likely once 
(1 — e)N/c cells have failed, since the initially failed cells reach their dynamical failure 

i 1/2 



1/2 

stress. The mean event size is roughly equal to n cf , therefore the mean time between 



events is TqtiJj /N, where To = (r Sj j — r a ^)/(KLv) is the characteristic time over which a cell 
is loaded from its arrest stress to its failure stress. The typical waiting time to see a switch 
to the runaway phase yields [cf. Eq.fli~5Dl 

T f ^T^ ^-w-;+* »\, («,) 

n c ' f { c n cf J 

with Cf r a factor of order unity which varies weakly with e and c in the region of interest, 
provided that each attempt is statistically independent of each other. We verified that in 
our simulations indeed no time correlations of event sizes were present, out to many times 
the characteristic time To (see also the discussion section). The distribution of persistence 
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times should then obey Poisson statistics with mean Tf_> r . Figure f| depicts the distribution 
of persistence times (with a fit to Poisson statistics) for N = 400 and c = 0.73. Mean 
persistence times depend very sensitively on the conservation parameter c, as shown in the 
inset of Fig. | @. 



B. Runaway Phase 

Let us next consider the runaway phase. Immediately preceding the first cell failure 
after a runaway event, the stress gaps {Ss n = — X n } have the probability distribution 
Prob(5s n — s) — exp [—pNs/{\ — e)] . Hence, {X n } has a mean fix(n) = n(l — e)/(pN) 
and variance cr\{n) = n(l — e) 2 /(pN) 2 . As long as dynamical effects are absent, the stress 
redistributed to each cell following the failure of the first n cells is still given by {Y n } with 
mean /iy(n) = nc/(pN) and variance 0y(n) w n(cW/pN) 2 <C o 2 x (n). Thus, the mean 
and variance of the step size for {Z n } is p — (1 — e — c)/(pN) and a 2 « (1 — e) 2 /(pN) 2 . 
The probability for an earthquake to terminate after n cell failures is [including finite size 
corrections in analogy with Eq. ([Op] 

r \ . A r f n(l + n/iV) ] . . 

p M (n)K _ ; _ exp |_i_LJ|, (17) 

2 



2(l-e 



b-(i- e )] 2 



{l + 0([c-(l- e )] 2 )}, Ol-e. (18) 



Since p, < for 1 — e < c, Z n < with finite probability for all n and a runaway event occurs. 
In fact, a runaway event is inevitable since Z^- < and the runaway event will commence 
once {Z n } reaches its maximum. The total number of cells that fail before a runaway event is 
given by the position of the maximum of {Z n }, whose probability distribution is proportional 
to n cr p^\n) for n 3> n cr . The mean number of these "precursor" cells is of order nj/ 2 , 
which remains a finite constant as iV — > oo, i.e., for big systems almost all the slip happens 
during the runaway events [|16j. 

The remaining cells will all fail during the runaway event. Imagine a situation where a 
fraction r > (1 — e)/c of the cells have failed. At that point, the total redistributed stress 
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per cell, is 



S = c 



r+(r- - [c + c 2 + ■■•} 



c 

c[r-(l-e)] 
1 -c 



(19) 



where the second term arises from repeat failures of some cells. S > 1 is needed to ensure 
that small event cells fail again and recreate the stress distribution (|H]). This is achieved if 

^ * 1 
r > r = e. 

c 

Thus, the large event cannot recreate the stress distribution fllPf ) if more than (1 — r*)N = 
r c N cells fail during background activity. This usually leads to a breakup of the bunched 
stress distribution and subsequent evolution towards the G-R stress distribution (§). The 
typical persistence time of the Runaway phase before a switch to the G-R phase is Jl5| 



T „, „ To ^l exp ( < C - + < C - WW] , e>C, (20) 
n 2 cr { c*cn cr J 

provided that all attempts are statistically independent of each other. (We have explicitly 
checked in the simulations that in the runaway phase, the particular realizations of stress 
distributions immediately following a large event are statistically independent of each other 
PI). T r ^f becomes comparable to the typical time between runaway events when c \ c*, 
as expected [0 . Figure § depicts the distribution of persistence times and a fit to Poisson 



statistics for N = 100 and c = 0.73. The inset shows the dependence of mean persistence 
times on c for iV = 100. Although agreement with Eqs. ([l6|) and (p0| ) is rather poor, the 
strong exponential dependence as a function of conservation parameter is evident. 



IV. DISCUSSION 

The persistence times in both the G-R phase and the runaway phase diverge exponen- 
tially with system size for (1 + e) -1 < c < 1, and the system remains in either phase for 
extremely long times, thus the phase space has two almost stable attractors. Clearly, the 
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runaway phase represents a more "ordered" stress distribution. Indeed, the basin of attrac- 
tion for the runaway phase is extremely small. In order to quantify this aspect, consider the 
"configurational entropy" for a given stress distribution p(s), with Sj = (rj— T a> i) / (rfj— T a ,i) ■ 



For the G-R phase, = 0, indicating that a "generic" stress distribution characterizes 

the G-R phase. On the other hand, in the runaway phase 



indicating that the stress distribution is highly organized in that phase. For discrete N the 
stress distribution is approximated by a histogram of the stress values, and the integral is 
replaced by the sum over all bins of the histogram. 

The time evolution of the configurational entropy of the stress distribution, calculated 
with a 10 bin histogram, is depicted in Fig. ^| along with event sizes. It is clear that S con f{t) 
can be used as an "order parameter" , a number that distinguishes the G-R phase and the 
runaway phase, with the advantage that it can be determined at any instant. Histograms of 
event sizes require a finite time interval to collect, and there is always the danger of mixing 
events from one phase with the other, thereby confusing the picture: The cumulative event 
size distribution over many persistence times is a weighted average of two entirely different 
event distributions, which obscures the underlying physical phenomena. S con f provides a 
reliable way to separate the two phases and makes it possible to accumulate accurate event 
size distributions for both of them. Unfortunately, such a quantity cannot be determined 
from existing field data since the spatial distribution of stress is unknown. 

So far, the discussion has centered around the If <C 1 limit, and the main role played 
by the heterogeneities has been the "randomization" of the stress distribution at time scales 
over which all cells fail a few times. The distribution of loading times, over which the 
cells are loaded from the individual arrest stress to the failure stress, has a mean T = 
{t s ,% — T a,i) I '{Kiv) and standard deviation of order WTq. Therefore, the "randomization" 




(21) 




(22) 
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time, over which the stress variables Sj become roughly uncorrelated, is of order Tq/W. 
Thus, even for small W, for large enough N this will be small compared to the persistence 
times, which scale exponentially in N [See Eqs. (|i~6D and (pUP]. This ensures the consistency 
of the assumption that the earthquakes are basically statistically independent of each other. 
The validity of this assumption of statistical independence can be explicitly verified by 
examining the time correlations of event sizes numerically; indeed no trace of any correlation 
was found in our simulations, out to many times the randomization time Tq/W. Likewise, 
we have explicitly checked that in the runaway phase, the particular realizations of stress 
distributions immediately following a large event are statistically independent of each other 

For finite values of W < 2, we expect most of the features to remain qualitatively 
unchanged: In the G-R phase, the exponential cutoff size still diverges as n c f ~ (1 — c)~ 2 , and 
although c* in general depends on W and the shape of p{s a ), there is still a persistent runaway 
phase for c* < c < 1. However, the situation is likely to change qualitatively once arrest 
stresses can be arbitrarily close to failure stresses, i.e. W = 2, and new values for the arrest 
stresses are picked every time a cell fails |TjJ. This corresponds to the situation discussed 
in Ref. 0] for finite-dimensional systems. Immediately upon introduction of dynamical 
weakening (e > 0), n c / ~ e~ 2 when c approaches 1, i.e. the cutoff size no longer diverges. 
Furthermore, for c = 1 the G-R phase is no longer persistent since the persistence time Ty_> r 
remains finite for large N. 

Some of the results presented for the mean-field model, especially the qualitative phase di- 
agram, calculated exponents for the power law earthquake distributions, and the divergence 
of the cutoff length scale, can be expected to apply to models with realistic interactions, up 
to logarithmic corrections. Such is because the underlying critical points that control these 
exponents remain mean-field-like down to 2 dimensional faults. This result is firmly estab- 
lished for the e = case @>HI- At finite e one expects the nucleation size for the runaway 
phase, which equals (1 — e)N/c in mean field theory, to become independent of the system 
size, since elastic forces in the fault plane concentrate stresses along the earthquake rupture 
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front as the earthquake progresses. Earthquakes bigger than a finite nucleation size N crac k 
become unstoppable in the presence of dynamic weakening effects and small disorder [|TJ, 
and rupture the entire fault. Nevertheless, for n c f < N crac k, the mean-field scaling results 
may still apply at finite e, provided that W < 2, i.e., there is a finite minimum stress drop 
associated with each cell failure. For systems with N crac k > N, this range will extend all 
the way to the fault size. In this case, one remarkable consequence is that since generically 
(1 — c) ~ 1/y/N HH, the cutoff size in the G-R phase n c / ~ (1 — c)~ 2 ~ N, i.e., earthquakes 
on individual fault zones obey power law statistics for events up to a finite fraction of the 
entire system size. 

An important result is the possibility that a fault system might switch from a 
"Gutenberg- Richter" earthquake distribution spontaneously to a "characteristic" earthquake 
distribution, as in the mean-field model. We note that calculations based on an entirely 
different model, simulating the coupled evolution of regional earthquakes and faults in a 
Theologically layered 3D solid ||12| , show similar behavior. Clear observation of such mode 
switching in nature requires data sets spanning many thousands of years. Paleoseismic stud- 
ies attempt to construct long histories of seismic events at given locations from sequences of 
displaces and highly disturbed rock layers. Remarkably, the longest available paleoseismic 
records, documenting large earthquake activity along the Dead sea transform in Israel || 
appear to be characterized by alternating phases of intense seismic activity lasting a few 
thousands of years, and periods of comparable length without large seismic events. Other, 
qualitatively similar alternating deformation phases have been documented in the eastern 
CA shear zone |Hj and the Great Basin Province in the western US JTl| . 



Another intriguing possibility might arise in a fault system of weakly coupled segments 
driven under similar conditions. The seismic response of such systems might exhibit a sort of 
"coexistence", i.e., a fraction of the patches might follow characteristic scaling whereas the 
others obey Gutenberg-Richter scaling, giving rise to a hybrid event size distribution. This 
may explain examples in the data of reference ||, where the characteristic "bump" in the 
distribution was not very pronounced. Finally, we note that part or all of the low magnitude 
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seismicity in the G-R phase may be too small to be detected by a seismic network. In this 
case the spontaneous switching between the Runaway and G-R phases may be interpreted 
as transitions from seismic response of a fault system to creep-like behavior. 
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19 



FIGURES 

FIG. 1. Schematic phase diagram of the system. There is a "coexistence" of two persistent 
stationary states called Gutenberg-Richter and Runaway phases, in a finite region of parameter 
space, marked region "(2) meta-stable" . For region 1 given by c < c* = 1/(1 + e) (line AB) one 
finds only small avalanches, i.e. the system is always in the Gutenberg-Richter phase. 

FIG. 2. Histograms of event size distributions in the two stationary states (phases), for 
W = 2/19, e = 0.5, N = 400. (a) The "Gutenberg-Richter" phase, characterized by a power-law 



earthquake distribution with an exponential cutoff. Solid lines are fits to the analytic form (15) 



with n c f as a fitting parameter. Also indicated are analytic estimates n£f ■ The inset shows a 
typical stress distribution of this phase for c = 0.7. The solid line is a fit to the analytic form (||). 
The nonuniform region near s = extends from —W/2 to W/2. (b) The "runaway" phase, with 
a similar background distribution and large characteristic events. The inset shows a typical stress 
distribution for c = 0.8. The solid line is a fit to the analytic form fllpp. The nonuniform region 
near s = extends from — W/2 to W/2. Near s = e it extends from e(l — W/2) to e(l + W/2). 

FIG. 3. The process {Z n }, which shows the incremental amount of stress needed to keep an 
earthquake going (see text for the precise definition). Each failure event corresponds to a segment 
of the process that starts out from a maximum up to that point and ends when it exceeds that level, 
and is marked as alternating circles and squares. The sample shown here, which corresponds to the 
stress distribution shown in the inset of Figure 2(a), depicts events of size 6, 14, 2, 1, 1, • • • . The fault 
is loaded adiabatically between these events, during the intervals when {Z n } moves monotonically 
up from one maximum to the next. These are shown as dotted lines connecting consecutive events. 
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FIG. 4. Distribution of persistence times T/_ r and T r ^ f for W = 2/19, e = 0.5, N = 100, 
c = 0.73. The lines are fits of the cumulative probabilities to the Poisson distribution. (Simulations 
for systems with other parameters that allowed for many more switches during the simulated 
times clearly also gave Poisson distributions for the distribution of persistence times.) Inset: The 
dependence of persistence times on conservation parameter c, (triangles Tf-> r , circles T r _>/) for the 
same values of W, e, and N. Statistical errors are comparable to symbol sizes. 

FIG. 5. Sample time series of earthquake sizes (top), plotted together with the conformational 
entropy S conf (t) (bottom) for W = 2/19, e = 0.5, N = 100, c = 0.73. (For the calculation of 
S con f(t), the simulated stress distribution was approximated by a 10 bin histogram of the stress 
values, which was evaluated immediately after each earthquake.) The earthquake size distribution 
changes drastically every time S con f toggles from to ln(l — e), indicating a transition from one 
phase to the other. A failed switching attempt from the G-R phase to the runaway phase is seen 
at about T/T = 10900. 
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